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A discrete Hubbard-Stratonovich transformation is presented for systems with an orbital degen- 
eracy A'^ and a Hubbard Coulomb interaction without multiplet effects. An exact transformation is 
obtained by introducing an external field which takes A'' -I- 1 values. Alternative approximate trans- 
formations are presented, where the field takes fewer values, for instance two values corresponding 
to an Ising spin. 
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Systems with orbital degeneracy have recently at- 
tracted much interest. For instance, in systems with 
colossal magnetoresistance, like the manganites, the de- 
generate 3c? orbital plays a major role. Another example 
is the alkali-doped Ceo compounds, A3 Ceo (A= K, Rb), 
where the partly filled ti„ orbital has a three-fold degen- 
eracy. Other well-known examples are the 3d metals. 

Quantum Monte Carlo, auxiliary field methods are 
popular methods for treating strongly correlated systems. 
In this approach a Hubbard-Stratonovich transformation 
is used to convert the many-body problem into a one- 
body problem—at the cost of introducing a fluctuating 
auxiliary field.EJ A substantial simplification was intro- 
duced by Hirsch, who showed that for a system with- 
out orbital degeneracy the field only needs to take tw| 
values and that it can be described by an Ising spin 
This method can be generalized to a system with the or- 
bital degeneracy TV by introducing an Ising spin for each 
pair of orbitals. This leads, however, to a large number 
{2N {2N —1) /2) of Ising spins for each site. This approach 
is very general, and can, for instance, handle the case 
when multiplet effects are considered. In many cases, 
however, we are interested in models which only have a 
Hubbard Coulomb interaction U , without any multiplet 
effects. This simplifies the problem substantially, since 
the Coulomb energy then only depends on the total num- 
ber of electrons on a given site, and not on the precise 
occupancy of the different levels. It is then natural to 
focus on the occupation number operator n for a given 
site. To calculate the partition function, we would then 
like to find an expression of the type 
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where n = 0, 1, 27V corresponding to the possible oc- 
cupancies of a given site. At is the step length in r 
introduced in the Trotter procedure and = exp(— x/j). 
Because there are 2iV-|-l conditions, it should be possible 
to satisfy these conditions exactly by letting k run over 
-I- I values, since there are then 2N + 2 parameters. 



Below we find that this is indeed the case. For U < 0, 
i.e. an attractive interaction, Wk is real and positive and 
Xk is real. For U > the Xfe's are complex, but come in 
pairs of complex conjugates. 

This problem can be solved analytically by rewriting 
the identity used by HubbardEl 
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where the weight function ^/{z) is defined appropriately 
and u — —UAt. Since the integral (^) now contains a 
polynomial of a maximum power 2N, the integral can 
be replaced by a discrete sum of the type in Eq. (|l|) by 
using the theory of Gauf5ian integration and orthogonal 
polynomials .B The transformation (|l|) is then exact for 
< n < 2A^ -I- 1. Below we reformulate this procedure 
in such a way that it can be generalized to a more useful 
form. 

We introduce a polynomial 
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where uq = 1- We also define 

B{n) = e-t^Ar«(„-i)/2^ 



(4) 
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which is the left hand side of Eq. (|^) . From Eq. (||) to- 
gether with P{zk) = we then obtain 
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where K = N +lior the moment, An^m = B{N '\-n — m) 
and Cn — —B{N + n). The solution of Eq. (H) gives a™. 
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These coefficients define P{z) from which we obtain the 
roots Zfe. The first iV + 1 equations in Eq. (|l|) then form 
a set of linear equations in Wk, which is solved. This 
approach exactly reproduces the result of the Gaufiian 
method mentioned above. Finally we transform from Zk 
to Xk by using Xk — — ln(2:fc). Since B{n) is real, the coef- 
ficients am are also real. The roots Zk of the polynomial 
P{z) are then real or come in pairs of complex conju- 
gates. The same is true for Xk- For U < Wk is positive 
and Xk is real, while for U > Wk and Xk in general 
come in pairs of complex conjugates. In the following we 
consider a positive (repulsive) U. 

By using K = N +1 in Eq. we really impose a phys- 
ically irrelevant condition, since it means that Eq. (Q) is 
satisfied for the occupancy n = 2N + 1 , which never oc- 
curs in the problem. We therefore relax this condition 
and use K = N, thereby satisfying Eq. (|l|) exactly for 
n = 0, ...,27V. Thus we put Atv+i,™ = i^TV+i,™- The 
value of Cn+1 = ctN+i can then be used to impose some 
additional condition on the transformation in Eq. (|^), as 
is discussed below. 

To use these results in a Monte Carlo approach, we 
introduce 



(8) 

where Pk > and J2kPk — 1- Wk are real and positive 
yk = 1, since J2j = 1- This leads to 
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FIG. 1. cr(n) as a function of Qjv+i for = 3 and 
f/Ar = 0.05. Results are shown for n = 3 (solid), n = 2,4 
(dotted), n = 1,5 (dashed) and n = 0, 6 (dashed-dotted). 



TABLE I. Values of Wk and xu (Eqs. (0,|l|)) and a{n) 
(Eq. (0)) for iV = 1 and no = iV as a function of Arf/. 
1^2 ~ wl and X2 — xl. The first lines were obtained by using 
K = N + 1 and the following results by using K — N and 
minimizing a{l){= 0). 



AtU 


Wl 


Xl 


a(0) 


a(l) 


a{2) 


0.1 


.5+i.0781 


.10-i.3097 


0.501 


0.156 


0.156 


0.2 


.5-Hi.l089 


.20-i.4290 


0.701 


0.212 


0.212 


0.4 


.5+i.l499 


.40-i.5823 


1.194 


0.300 


0.300 


0.1 


.5 


.05-1.31360 


0.324 


0.000 


0.324 


0.2 


.5 


.1 -i.43980 


0.471 


0.000 


0.471 


0.4 


.5 


.2 -i.61160 


0.701 


0.000 


0.701 



We now use a probability interpretation and choose the 
term yfee"^*"" with the probability pfc. We introduce the 
relative standard deviation 

t-[ ^(") 

The free parameter a^+i can now be used to minimize 
a{n). For instance, for a half-filled systems with a large 
value of U, most sites have the occupancy N. It should 
then be useful to minimize cf{N). Fig. |l| shows cr(n) as 
a function of a^+i for = 3. It illustrates that it is 
possible to obtain a{n = 3) = and at the same time to 
also obtain rather small values for n — 2 and n = 4. If 
the system instead is close to some other integer filling 
n, we can choose ciAr+i so that the corresponding cr(n) is 
small. 

Since we are often interested in a system at or close 
to some integer occupancy uq, it is also convenient to 
introduce 



where Wk is related to Wk and Xk by a simple transfor- 
mation 

Table | shows results for Wk, Xk and a{n) for N = 2 
and hq = N. The first lines were obtained by satisfying 
Eq. (0) also for the irrelevant case n = 2 + 1 , while the 
following lines were obtained by satisfying Eq. (^ only 
for n < 2N and using the additional freedom to mini- 
mize <j{N). In the latter case cr(n) is symmetric around 
n = N and Wk (but not Wk) is real. 

For = 1 and = 2 we have obtained analytical ex- 
pressions for the coefficients lik and Xk which minimize 
(j{N). Using no = N we find the exact transformation 
for A^ = 1 

wi^2 = 0.5; xi^2 = 0.5UAt ± itan" Ve*^"^^ - 1. (13) 

This result is identical to a result obtained by IIirsch,i 
but presented in a different form. For A^ = 2 we obtain 
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wi,2=a; a;i,2 = 1.5i7AT±icos"^(2/3- 1) 

W3 = l- 2a; X3 ^ 1.5UAt (14) 



l)/4, 7 = exp(-0.5C/AT) and 
We observe that the w's are all 



with /3 = - l)/(7 
a = (7 - l)/(2/3 - 2)/2 
real in these cases. 

For N ~ 3 we have not found a simple analytical ex- 
pression for lik and Xk- Therefore numerical results are 
eriven m Table || for no — N. The corresponding values 



of (T(n) are given in Table III 



TABLE II. Values of Wk and Xk for A'" = 3 as a function of 
AtU obtained by minimizing o-{N). {x2=x'l and X4=x'^). 



AtU 




W3 


Xl 


X3 


0.05 


.45154 


.04846 


.1250-1.16385 


.125 -1.51551 


0.1 


.44890 


.05110 


.25-i.22883 


.25 -1.72002 


0.2 


.44345 


.05655 


.5-1.31561 


.5-1.99334 


0.4 


.43204 


.06796 


1.-1.42465 


1.-11.33776 



TABLE III. Values of a(n) for A = 3 as a function of AtU 
obtained by minimizing cr{N) and corresponding to the Wk's 



AtU 


a(0) 


a(l) 


a(2) 


a(3) 


a(4) 


a(5) 


a(6) 


0.05 


0.75 


0.47 


0.23 


.0000 


0.23 


0.47 


0.75 


0.1 


1.21 


0.70 


0.32 


.0000 


0.32 


0.70 


1.21 


0.2 


2.25 


1.11 


0.47 


.0000 


0.47 


1.11 


2.25 


0.4 


5.97 


1.99 


0.70 


.0000 


0.70 


1.99 


5.97 



So far we have only considered exact transformations 
of Eq. (|l|), which requires a field with N+1 values. Since 
the sampling procedure as well as the Trotter decomposi- 
tion introduce errors, we may also consider approximate 
transformations. This has the advantage that we can 
then reduce the number of values of the field and make 
the sampling procedure more efficient. For this purpose 
we introduce the relative error 
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should be minimized, where a is some appropriate pa- 
rameter which emphasises the suppression of the errors 
A (a << 1) and f{n) is a weight factor which empha- 
sizes the terms close to n = A^. The minimum of Eq. ( p^ ) 
is searched using a Metropolis algorithm followed by a 
Newton minimization. The precise values of Wk and Xk 
depend on the weight factors a and /,„ . There is also no 
guarantee that the Metropolis search finds the absolute 
minimum. 



In Table iV we show results for A'' = 2 using a field 
with just two values. The table shows that if UAt is 
not too large, the errors A are small. For instance, for 
UAt = 0.4 the largest error is just 0.04, which occurs 
for n — 1 and ri = 3. Such an error means that we are 
effectively using a Hamiltonian where U depends on n 
and has an error of about 4 % for n — 1 and n — 3, but 
is almost exact for other values of n. It should therefore 
be possible to treat N = 2 with just a single Ising spin, 
as for A^ = 1, although the field in the A'' = 2 is complex 
in contrast to tbp real field in one of the Hirsch A^ = 1 
transformations. a Table ^ shows results for A^ = 3 and 
just two values of the field. In this case the errors are 
larger than in the N — 2 case, but the larger errors hap- 
pen for the configurations with n = and n = 6, which 
should be rare for large values of U. 

TABLE IV. Values oi w, x, A and cr for A = 2 as a func- 
tion of AtU obtained by tfie expression in Eq. ( p^ using just 
two values of the field. A(n) and (T(n) are symmetric around 
n = 2. 



AtU 


Wi 


Xl 


A(0) 


A(l) 


A(2) 


a(0) 


a(l) 


a(2) 


.05 


.5 


.075+1.21990 


.0000 


.0006 


.0000 


.47 


.22 


0.0 


.1 


.5 


.15+1.30580 


.0000 


.0025 


.0000 


.70 


.32 


0.0 


.2 


.5 


.3+1.41808 


.0000 


.0100 


.0000 


1.11 


.45 


0.0 


.4 


.5 


.6+1.55239 


.0000 


.0397 


.0000 


1.99 


.64 


0.0 



TABLE V. Values of A{i) for A = 3 as a function of AtU 
obtained by the expression in Eq. ( p^ using a field with just 
two values of the field. A(n) and (T(n) are symmetric around 
n = 3. The values of (T(n) are similar as in Table 



III 



AtU 


Wl 


Xl 


A(0) 


A(l) 


A(2) 


A(3) 


0.05 


.5 


.125+1.21990 


-.0104 


.0000 


.0006 


.0000 


0.1 


.5 


.25+1.30580 


-.0466 


.0000 


.0025 


.0000 



In the same spirit we can require that Wk and Xk are 
all real. It is still possible to obtain parameters so that 
Eq. (|l]) is rather well satisfied. For an equal number of 
spin up and spin down electrons, even away from half- 
filling, the spin up and spin determinants are then equal 
and real and their product is positive definite. Thus the 
determinants do not cause a sign problem in this case. 
For a repulsive U some of the Wk must, however, be neg- 
ative. The sign problem then enters in a different form. 
It would be interesting to test whether or not this sign 
problem is less serious than in other formulations, where 
it enters via the determinants. 

To summarize, we have presented an exact discrete 
Hubbard-Stratonovich transformation for a system with 
the orbital degeneracy A^ using a field which takes A^ + 1 
complex values. The sign problem in the traditional 
Quantum Monte Carlo treatment is therefore converted 
into a phase problem. The standard deviation of the 
transformation has been minimized. It will be interesting 
to see how this influences the sampling efficiency. Ana- 
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lytical formulas for the fields were provided for iV = 1 
and N = 2 and numerical results for TV = 3. Alternative 
approximate transformations were presented, where the 
number of values of the field is smaller than N + 1. For 
instance, it was found that for A'' = 2 and N = 3 and 
for AtU not too large, a rather accurate transformation 
can be obtained with a field which only takes two values, 
corresponding to an Ising spin. 

We would like to thank A. Muramatsu for very stimu- 
lating and helpful discussions. 
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